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Abstract 

The relative stability and ergodicity of deterministic time-reversible thermostats, both singly and 
in coupled pairs, are assessed through their Lyapunov spectra. Five types of thermostat are coupled 
to one another through a single Hooke’s-Law harmonic spring. The resulting dynamics shows 
that three specific thermostat types, Hoover-Holian, Ju-Bulgac, and Martyna-Klein-Tuckerman, 
have very similar Lyapunov spectra in their equilibrium four-dimensional phase spaces and when 
coupled in equilibrium or nonequilibrium pairs. All three of these oscillator-based thermostats 
are shown to be ergodic, with smooth analytic Gaussian distributions in their extended phase 
spaces ( coordinate, momentum, and two control variables ). Evidently these three ergodic and 
time-reversible thermostat types are particularly useful as statistical-mechanical thermometers and 
thermostats. Each of them generates Gibbs’ universal canonical distribution internally as well as 
for systems to which they are coupled. Thus they obey the Zeroth Law of Thermodynamics, as a 
good heat bath should. They also provide dissipative heat flow with relatively small nonlinearity 
when two or more such bath temperatures interact and provide useful deterministic replacements 
for the stochastic Langevin equation. 
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I. DETERMINISTIC TIME-REVERSIBLE THERMOSTATS AND ERGODICITY 


In the early days of equilibrium many-body simulation, intercomparisons of results from 
constant-temperature Monte Carlo with those from constant-energy molecular dynamics 
were indirect, requiring the conversion of isothermal NVT data to isoenergetic NVE data. 
Only for hard disks and hard spheres, where temperature and energy are proportional, are 
the two ensembles identical^. 

Shuichi Nose provided a conceptual breakthrough linking temperature T and energy E 
by discovering an isothermal canonical-ensemble molecular dynamics^i^. He started out by 
including an additional time-reversible “time-scaling variable” s in his novel Hamiltonian : 

T^Nose = [ K{p)/S^ ] + ^q) + {pI/2M) + (#A;T/2) ln(s2) . 

Here is the number of degrees of freedom, including s , and ps is the momentum conjugate 
to s . M is a free parameter. Getting to the canonical ensemble from Nose’s Hamiltonian 
equations of motion requires just two more steps ; [ i ] “scaling the time”, multiplying all of 
Nose’s time derivatives by s : 

{ q ^ sq p ^ sp} ] s ^ ss ] Ps ^ sps ] 

then [ ii ] replacing all the scaled momenta { {p/s) } by { p } . Nose showed that the 
resulting distribution for the { g,p } is Gibbs’ canonical distribution. 

Starting instead with the “Nose-Hoover” equations of motion4^- , 

#=ND 

{q = {p/m) ; p = E{q) - (p } ; ( = ^ [ {p^/mkT) - 1 ]/r^ , 

Hoover showed that Gibbs’ canonical distribution, with s absent and with ^ = ND not 
including that extraneous s variable, is a stationary solution of the extended phase-space 
continuity equation , 

it it 

{df/dt) = 0 = - Y.{dfq/dq) - Y.{dfp/dp) - {dfC/dQ . 

Here C is a “friction coefficient” proportional to Nose’s Ps and the free parameter r takes 
the place of Nose’s parameter M . The stationary distribution function for the NosGHoover 
motion equations is canonical in the { O', p } and Gaussian in the friction coefficient C, : 

f{(l,P,C) oc ; K = ^{q) + K{p) . 
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The relaxation time r controls the decay rate of velocity fluctuations in a D-dimensional 
A^-body Hamiltonian system. 

In order for averages using the Nose-Hoover equations to agree with canonical-ensemble 
averages, it is necessary in principle that the dynamics be “ergodic”, meaning that it must 
reach all of the { q,p } phase-space states. Of course in practice only a representative 
sampling of such states can be achieved. Whether or not the motion equations are “ergodic” 
in this sense, capable of reaching all of the states described by Gibbs’ canonical ensemble, 
depends upon the details of the underlying Hamiltonian 'H{q,p) . For a simple harmonic 
oscillator Hoover pointed out that neither the four original unsealed Nose equations nor the 
three NosGHoover motion equations are ergodic : 

q = +(p/'5^) ; P = —q ] Ms = Ps = [ /s^) — (2/s) ] ; [ Nose, Not Ergodic ] . 

q = +P ; P = “ Cp j C = [ “ 1 1/'^^ ; [ Nose — Hoover, Not Ergodic ] . 

These Nose equations are singular while the Nose-Hoover equations are not. 

With the relaxation time r equal to unity, numerical NosGHoover calculations show 
that about hve percent of the initial conditions drawn from the Gaussian distribution 
g-'?V2g-p’^/2g-cV2 chaotic, making up a “chaotic sea” in the {q,p,C) space which is 
penetrated by an inhnite number of holes. The remaining 95% of initial conditions gener¬ 
ate regular periodic toroidal solutions which £11 in these holes. Evidently this NosGHoover 
thermostated oscillator is far from “ergodic” [ where in this case an ergodic solution would 
have a smooth analytic Gaussian density throughout {q,p,C) space ] : 

/Gibbs(<?,P,C) oc [ Ergodic ] . 

In 1990 Bulgac and Kusnezov reiterated the usefulness of the phase-space continuity 
equation in formulating more complicated motion equations ( with two or three thermostat- 
ing control variables ). Shortly thereafter they were able to use this approach to fill out 
the complete Gibbs’ distribution for prototypical Hamiltonian systems like the harmonic 
oscillator and the two-well “Mexican Hat” problemsS.^— . 

There is no foolproof test for ergodicity. The only reliable diagnostic for space-filling 
ergodicity is the Lyapunov spectrum^^— . If this spectrum, which measures the long-time- 
averaged sensitivity of the dynamics to initial conditions, is not only chaotic, but also the 
same for all initial conditions, the system is likely ergodic. For the harmonic oscillator 
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this means that all {q,p) oscillator states, all the way to infinity, will eventually occur. 
In nonergodic systems the long-time-averaged spectrum depends upon the initial condi¬ 
tions. In ergodic systems the long-time spectrum is always the same. By 1996 there were 
three known simple types of differential equations [ HH = Hoover-Holian, JB = Ju-Bulgac, 
MKT = Martyna-Klein-Tuckermann ] which were thought to produce ergodicity in a one¬ 
dimensional harmonic oscillator-'i^“— . In addition to propagating the phase variables (g,p) 
all three included two additional control variables (C^O capable of generating Gibbs’ canon¬ 
ical distribution. For the simple harmonic oscillator these three models take the form : 

{q=p] p =-q-Cp-ip^ ] C = p^ -T ; i = p'^ - 3Tp^ } [ HH ] ; 

{q=p; p = -q - ; ( = p‘^ - T ; ^ - STp^ } [ JB ] ; 

{q = p- p=-q-(p- ( = p^-T-^C; ^ = C' - T } [ MKT ] . 

So far no one-thermostat oscillator system ( with just a single control variable ) has been 
shown to be ergodic though it certainly is possible that such a one will be discovered. 
We explicitly include the temperature T in all of these models in order to set the stage 
for nonequilibrium systems incorporating both “cold” and “hot” degrees of freedom. It 
is worth pointing out that not all two-thermostat oscillator systems are ergodic. Patra 
and Bhattacharya showed that their very reasonable set of doubly-thermostated oscillator 
equations is not ergodio^^ii^ : 

{q = p-^q- p=-q-C^p- C^=p^ -T ■ ^ = g2_T}[PB]. 

Our work here has three different aspects. First we explore the equilibrium Lyapunov 
instability which facilitates the ergodicity of all three thermostats. In the equilibrium case, 
where the phase-space distribution is a smooth Gaussian, the Lyapunov instability embedded 
in that Gaussian has a totally different multifractal character. We show that the equilibrium 
HH, JB, and MKT thermostats obey the Zeroth Law when coupled with one another or with 
other Hamiltonian systems. These thermostats can provide Gibbs’ complete canonical dis¬ 
tribution provided that any internal energy barriers are not too large relative to kT. Second 
we consider nonequilibrium cases, where heat flows between thermostats with the thermostat 
temperatures set at different values. The nonequilibrium dynamics is still ergodic but the 
phase-space distribution is no longer a smooth Gaussian. It is instead an intimate multi¬ 
fractal combination of the attractor-repellor pairs common to macroscopic time-reversible. 
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but dissipative, systems. We also illustrate the application of ergodic thermostats to kinetic 
barrier-crossing problems. Last of all, we summarize the lessons learned in this study of 
three different sorts of applications of our ergodic time-reversible thermostats. 


II. LYAPUNOV EXPONENTS-THERMOSTATED HARMONIC OSCILLATORS 

The most convincing evidence for the lack of ergodicity with one thermostat variable, 
and its lack or presence in two-thermostat systems, is the Lyapunov spectrum. Here we con¬ 
centrate on thermostated oscillators, the prototypical “difficult” case. A harmonic-oscillator 
system of four ordinary differential equations, such as the HH, JB, MKT, and PB sets for 
{QjpyCyO four such exponents (Ai, A 2 , A 3 , A 4 ) . These exponents describe the deforma¬ 
tion of an inhnitesimal hypersphere or hypercube in the four-dimensional phase space which 
contains the motion. Shimada and NakashimaA, as well as Benettin’s group>^ described 
general approaches to determining the spectrum of exponents. 

In an n-dimensional space their algorithms require the solution of n-|-l sets of n equations. 
In addition to a “reference trajectory” the equations describe the motion of an associated or¬ 
thonormal set of n comoving basis vectors, centered on the reference trajectory and locating 
n nearby “satellite trajectories” in the n-dimensional space. The hrst Lyapunov exponent 
gives the time-averaged rate at which two neighboring solutions of the equations diverge 
from one another, ( 5{t) ~ ) . The rate at which the area dehned by three nearby 

solutions ( the reference and two others ) diverges ( or converges ) ~ dehnes A 2 , 

while the rate at which the volume dehned by four solutions and the hypervolume dehned 
by all hve of them dehne A 3 and A 4 . Although these exponents are typically “paired” for 
Hamiltonian systems, with 

[ Ai = ( Ai(t) ) = —( \i{t) ) ; A 2 = ( A 2 (t) ) = —( A 3 (t) ) = 0 ] —t Aj = 0 , 

1 

with the sum total equal to zero, as implied by Liouville’s Theorem, none of the two- 
thermostat systems is Hamiltonian so that this instantaneous symmetry is missing in the 
time-dependent exponents { A(t) } . But because the equations of motion are time-reversible 
the time-averaged spectra { A } are symmetric. 

For an ergodic system the long-time-averaged spectrum of exponents should be indepen¬ 
dent of the initial conditions. In practice, for the harmonic-oscillator system, a few hundred 
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Figure 1: The sign of the local (time-dependent) value of the largest Lyapunov exponent for 
the Hoover-Holian oscillator is shown in the (g,p, 0,0) plane. Because a mirror reflection of the 
dynamics ( perpendicular to the q axis ) changes the signs of both q and p the Figure reveals an 
inversion symmetry. In the upper half plane red indicates a positive exponent and green a negative. 
The colors are reversed in the lower half plane. The probability density in the cross-sectional plane 
is a simple Gaussian in q and p . Nevertheless the local Lyapunov exponents vary in a multifractal 
manner throughout the four-dimensional space, reflecting their sensitivity to bifurcations in their 
past histories. In Figures 1-4 both q and p range from -3 to -|-3 . 

oscillations is time enough to indicate whether or not a chosen initial condition ( either 
from a grid or chosen randomly ) converges to a spectrum close to the unique long-time- 
averaged Lyapunov spectrum. Estimates of these long-time-averaged exponents, after a time 
t = 40 000 000 using a fourth-order Runge-Kutta integrator with an adaptive timestep, are 
as follows : 

{ }hh = -|-0.068o, -|- 0 . 000 , — 0 . 000 , —O.O 680 ; 
in the Hoover-Holian case, and 

{ A }jB = +0.0797, +0.0000, -0.0000, -0.0707 ; 

in the Ju-Bulgac case, and 

{ }mkt — +O.O 665 , +0.0000, —0.0000, —O.O 665 . 
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Figure 2: In the upper half plane red indicates a positive local Lyapunov exponent for the ergodic 
Ju-Bulgac oscillator model. The inversion symmetry ( with change of color ) is obeyed by all these 
time-reversible oscillator models. The probability density in the four-dimensional space is Gaussian 
in ( 7 , p, C, and as it is also in the (g,p, 0 , 0 ) plane displayed here. 

in the Martyna-Klein-Tuckerman case. 

The increase in the largest Lyapunov exponent with the number of quartic forces ( none 
for MKT, one for HH, and two for JB ) suggests, as emphasized by Bulgac and Kusnezov^, 
that these terms are particularly well-suited to promoting chaos and/or ergodicity. It seems 
likely that sextic forces, controlling ( p® ) , would have no particular advantages and would 
slow numerical work. 

The probability densities for the three cases follow easily from the phase-space continuity 
equation. This makes it is easy to check for ergodicity by comparing relatively short-time 
estimates for the Lyapunov spectrum starting out with initial conditions chosen according 
to the stationary multivariable Gaussian distributions. 

The ergodicity of the Martyna-Klein-Tuckerman oscillator was called into question by Pa- 
tra and Bhattacharya in 2014^^—, based on apparent “holes” in a (g,p, C, 0 = +1) 

double cross-section plane, the analog of a Poincare plane for a three-dimensional flow prob¬ 
lem. To resolve this question, the subject of the 2014 Ian Snook Prize^^, we chose one million 
different initial conditions from the appropriate four-dimensional Gaussian distributions for 
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the HH, JB, and MKT cases and determined that all of them were fully consistent with 
a unique chaotic spectrum. Each of these initial conditions was followed for one million 
fourth-order Runge-Kutta timesteps of 0.005 each. 

The alternative possibility, a regular ( nonchaotic ) solution with all four of the long-time- 
averaged Lyapunov exponents equal to zero, is easy to distinguish from the chaotic case. 
As an additional check the velocity moments generated by these three ergodic thermostat 
types likewise reproduce Gibbs’ values for the second, fourth, and sixth velocity moments 
with hve-figure accuracy : 

( pi 2,4,6 } ) = { 1 . 00000 , 3 . 0000 , 15.000 } . 

By contrast, the motion equations for a NosGHoover harmonic oscillator with a single 
friction coefficient ({t) and unit relaxation time , 

{ q = p ■ p = -q - (p ■ ( = p‘^ -I j ^ 

are known not to be ergodic despite their simple three-dimensional Gaussian solution 
of the phase-space continuity equation. Ghoosing initial conditions from this stationary 
distribution /nh oc and measuring the tendency of a nearby initial con¬ 

dition to separate gives the largest Lyapunov exponent for the chosen initial condition, 
Ai = ( {d/dt){ln{Sr)) ) . Because this system is conservative and time-reversible the long- 
time-averaged Lyapunov spectrum is symmetric, { -|-Ai, 0, —Ai } , and sums to zero. 10 000 
initial conditions followed for a time t = 50 000 divide neatly into two groups : 

0.00002 < Ai < 0.00014 ; 0.006 < Ai < 0.017 . 

557 initial conditions correspond to a Lyapunov exponent signihcantly different to zero while 
9443 correspond to nonchaotic toroidal solutions for which all three Lyapunov exponents 
vanish. Evidently roughly 6 % of the stationary measure is chaotic ( with a positive time- 
averaged largest Lyapunov exponent ). The remaining 94% consists of regular tori with a 
vanishing Lyapunov spectrum : 

( Ai, A 2 , A 3 ) = { 0, 0, 0 } . 

In the language of control theory the NosGHoover equations control the time-averaged 
value of the kinetic temperature, ( {p^/mk) ) = T . Subsequent work strongly suggests 




Figure 3: Red indicates a positive Lyapunov exponent in the upper half plane, negative in the 
lower. If the equations of motion are run backward ( by changing the righthandsides of all four 
equations ) the effect is to reflect the section shown here, changing the sign of q. Like the other 
models this Martyna-Klein-Tuckerman oscillator has a four-dimensional Gaussian distribution in 
its (g,p, C)0 phase space which the dynamics samples ergodically. 

that an ergodic oscillator requires control of at least two moments, not just one. Such an 
approach requires four or more ordinary differential equations. Bulgac and Kusnezov^ as 
well as Ju and Bulgao^ added a hfth equation so as to be able to thermostat a free particle 
( or a rotating and translating cluster of particles ) undergoing Brownian motion. 

The /onr-dimensional phase-space continuity equation shows that both the Martyna- 
Klein-Tuckerman and the Hoover-Holian thermostats are consistent with exactly the same 
four-dimensional Gaussian, /hh = /mkt • Unlike the Nose-Hoover single-thermostat model 
we believe that the three two-thermostat models are all ergodic. We believe that this 
proposition has been thoroughly conhrmed by the present work, conhrming the chaos of 
one million independent initial conditions in each case. In view of the lack of other suitable 
entries in the 2014 competition, we have awarded ourselves the Snook Prize for this hnding. 
It should be noted that all of these thermostat models, ergodic or not, are time-reversible, 
with the friction coefficients (C,0 the momentum p simply changing sign if the time- 
ordered sequence of coordinate values { q{t) } is reversed. For this reason the time-averaged 
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Figure 4: Red indicates a positive Lyapunov exponent in the upper half plane, negative in 
the lower. If the equations of motion are run backward by changing the righthandsides of all 
four equations the effect is to reflect the section shown here, changing the sign of q. Like the 
other models this Patra-Bhattacharya oscillator has a four-dimensional Gaussian distribution in 
its {q,p,Cj^) phase space which the dynamics samples ergodically. Unlike the other models about 
half the PB measure is regular rather than chaotic. 

Lyapunov spectrum is symmetric, with Ai -|- A 4 = A 2 + A 3 = A 2 = A 3 = 0 . Despite the 
simple analytic nature of the phase-space distribution, the chaos implied by a positive Ai 
engenders an amazing complexity, singular in its spatial behavior, to which we turn next. 

In exploring the HH, JB, and MKT ergodic oscillator solutions we uncovered an amazingly 
intricate multifractal structure most simply described in terms of the local largest Lyapunov 
exponent, Ai(t) . The three phase-space distributions { f{q,p, C, 0 } ^^^e all multidimensional 
Gaussians, of the form , 

where n is 2 for HH and MKT and 4 for JB. A closer look, using Lyapunov instability 
as a tool, reveals the interesting structures shown for the three thermostated oscillators in 
Figures 1-3 . For completeness, the corresponding Patra-Bhattacharya cross section is 
shown as Figure 4 . 

For each of the three ergodic models we show the sign of the local Lyapunov exponent. 
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Ai(t) at the location (g,p, 0,0) . The plots are analogous to a Poincare section for the 
more usual three-dimensional flows. At the equilibrium, T = 1 , these instability plots 
can be generated in either of two ways: [ i ] Whenever a (g,p, 0,0) trajectory comes close 
to (CjO = (0)0) plot the sign of the corresponding Lyapunov exponent at the location 
(g,p) ; [ ii ] Starting at (g,p, 0,0) integrate backwards in time for several hundred thousand 
timesteps, storing the backward trajectory. Then process the reversed trajectory data in 
the forward direction, with Lyapunov instability controlling a nearby constrained satellite 
trajectory, hnding the local Lyapunov exponent at the chosen endpoint location (g,p, 0,0). 
It turns out that these two methods are roughly comparable in cost. The second method 
has the advantage of providing more detail in regions where the measure is small. The 
precision shown in Figures 1-4 can be generated in simulations taking a day or two on a 
single processor. Method [ ii ] above, reversing a stored trajectory, is particularly well-suited 
to parallel simulations. The Figures use color to indicate the sign of the local exponent. 


III. FRACTALS LURK IN THE SIMPLE GAUSSIAN DISTRIBUTIONS 

These cross sections have a fractal look. In fact they are. The reason for this was 
pointed out in our paper with Dennis Isbister— . As the “past” history of a point is being 
generated, backward from the “present”, “bifurcations”, where the direction of the trajectory 
is uncertain at the level of the numerical work limit the trustworthy scale of the Lyapunov 
exponent when the stored integration is reversed, going forward in time. For the megapixel 
details shown in Figures 1-4 a few hundred thousand fourth-order ( or hfth-order or adaptive 
) Runge-Kutta timesteps suffice for reliable results at the resolution of the plotted sections. 
More detail can always be generated, at a cost exponential in the number of signihcant 
hgures, by reducing the numerical error of the trajectory integration. 

The cross sections shown in the Figures have inversion symmetry. As the source of this 
symmetry is not so obvious let us describe it. Consider a solution of the equations of motion 
going forward in time, { q(t),p{t),((t),^(t), Xi(t) } . Viewed simultaneously in a mirror 
perpendicular to the q axis an observer sees { —q(t),—p{t),(it),^{t),Xi{t) } . Evidently 
there is an inversion symmetry, with exactly the same Lyapunov exponent at two different 
values of {q,p) : 

Xi{+q,+p) = Ai(-g,-p) . 
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This observation saves a factor of two in compnter time. The two-dimensional cross sections 
where both of the friction coefficients vanish, C = = 0 are the starting points for the 

backward integrations. The Figures indicate the magnitude of the largest Lyapunov expo¬ 
nent obtained at the original starting point from a reversed trajectory, with red positive and 
green negative as the local exponents varies with position. This is the case for positive q . 
The colors are reversed for negative coordinate values in order to eliminate the discontinuity 
that would otherwise occur at vanishing q . 

As is usual the fluctuations in Xi{t) are an order of magnitude larger than its long-time- 
averaged value. The variation of the local Lyapunov exponent is not quite continuous. If the 
local exponent is plotted along a line in the {q,p) plane there are occasional jumps, present 
on all scales, in the exponent value. The fractal dimension accounts for the variation of the 
jump magnitudes as a function of resolution. In the MKT casein the data along such a line 
gave a fractal dimension of 1.69 rather than 1.00 . Studies of the jumps along such lines 
reveal their fractal structure, seen here for the first time in two dimensions. 

IV. THERMOSTATING CONFIGURATIONAL TEMPERATURE ? 

All three of the two-thermostat HH, JB, and MKT models for heat-bath control of a 
harmonic oscillator are ergodic. Because of this it was a complete surprise that Patra and 
Bhattacharya’s control of both the kinetic and the configurational^ oscillator temperatures, 
using two thermostat variables, was unsuccessful, leading to a mix of regular and chaotic 
solutions. Their model was : 

{ q = p - ■ p = -q - (p ■ ( = p"^ - I ■ ^ = g2 _ ^ | ^ ^ ) = (1, 1) [ PB ] . 

The stationary phase-space distribution function for the PB equations is the same as that 
characterizing the HH and MKT models : 

/hh = /mkt = /pB = . 

But unlike the HH and MKT equations the Patra-Bhattacharya model fails the Lyapunov- 
exponent test for ergodicity. About half the initial conditions chosen from the four¬ 
dimensional Gaussian distribution give regular rather than chaotic solutions. 

Although controlling coordinate moments like ( g^, g'^, g® ... ) may seem unphysical Lan¬ 
dau and Lifshitz showed ( in the 1951 Russian edition of their Statistical Physics ) that 
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what later came to be termed a “configurational temperature”—, based on a mean-squared 
force, can be derived from Gibbs’ canonical distribution. A single integration by parts, in 
the canonical ensemble, relates the mean-squared force to the curvature of the potential 
through the temperature T : 

kT j [ ]e-^/^^dq = J [ (VnY dq . 

So control of ( ) can be viewed as regulating the fluctuations in the forces. This concept is 

more appealing than coordinate control. But the lack of a physical “thermometer” capable 
of measuring force-based temperature is a disabling disadvantage of configurational temper¬ 
ature. The possibility of divergent or negative values of the instantaneous configurational 
temperature are further caution flags. The Campisi “thermostat” is an example of both 
defects : 

n = (pV 2) + (T/2) ln(g2) ^{q = p;p= -{T/q) } . 

Formally, Campisi dynamics is consistent with a Gaussian momentum distribution. Al¬ 
though at first glance appealing, this simple model exhibits a negative configurational tem¬ 
perature, —T , in one dimension and a divergent one in two^i. Disagreement between the 
kinetic and configurational temperatures is a clear violation of the Zeroth Law of Thermo¬ 
dynamics, to which we turn next. 


V. COUPLED THERMOSTATS-THE ZEROTH LAW OF THERMODYNAMICS 

To what extent do these minimalistic thermostats, described by one or two additional 
control variables, exhibit the characteristic properties of macroscopic thermal baths? An 
essential characteristic of such baths is their consistency with the Zeroth Law of Thermo¬ 
dynamics. Two “good” thermal baths, both at a temperature T, should, when coupled 
together or to another canonical equilibrium system at the same temperature T, remain at 
that common temperature without modifying Gibbs’ equilibrium distribution in either bath 
or in another equilibrium system. 

For our thermostated oscillators this desirable property can be tested by coupling pairs 
of heat baths together with a simple Hooke’s-Law spring. The modification of the equations 
of motion for oscillator number 1 and oscillator number 2, both at unit temperature, are as 
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follows : 


012 = (gi - ^ 2)^/2 -)■ 01 = 0i(T = 1 ) + g 2 - gi ; p 2 = P 2 {T = l) + qi- q 2 . 

It is easy to verify that the two Nose-Hoover oscillators, coupled together by a spring, 
have nonergodic solutions, corresponding to a six-dimensional analog of the tori making up 
most of the single-oscillator phase space. One such solution occurs with initial values of 
the two oscillator momenta {pi,P 2 ) = (0.99,1.01) . Accordingly we restrict our detailed 
investigation to the three ergodic two-thermostat heat baths to see whether or not they can 
obey the Zeroth Law. 

All the six combinations of like and unlike pairs of ergodic two-thermostat oscillators : 

[ HH+HH, HH+JB, HH+MKT, JB+JB, JB+MKT, MKT+MKT ] , 

provide chaotic eight-dimensional regions in phase space with symmetric Lyapunov spec¬ 
trum. This symmetry implies ergodicity. Figure 5 shows both the spectra and the spectral 
sums ( which give the growth rates of phase-space hyperspheres as functions of their dimen¬ 
sionality ). Even though there are relatively small but signihcant differences between the 
HH, JB, and MKT Lyapunov exponents pairing any two of these systems together ( with a 
Hooke’s-Law spring ) gives no indication of the dissipation which would result if one phase- 
space bath were to grow at the expense of the other. We conclude from these six examples 
that the HH, JB, and MKT heat baths all obey the Zeroth Law of Thermodynamics, quite 
a good outcome for dynamical systems which represent a thermodynamic heat bath with 
the low cost of only four phase-space dimensions. 

The largest Lyapunov exponent in these eight-dimensional problems varies from 0.20 for 
two coupled MKT oscillators to 0.23 for two coupled JB oscillators. The exponent for two HH 
oscillators lies in between. Evidently the quartic forces in the Ju-Bulgac and Hoover-Holian 
thermostats are slightly more effective in promoting chaos than are the cubic Martyna- 
Klein-Tuckerman forces. This same conclusion follows from the time required to generate 
the structures seen in Figures 1-4 . The additional chaos comes at the ( rather small ) price 
of an increased stiffness in the differential equations. If any of the three ergodic thermostats 
were not ergodic we would expect to see a dissipative strange attractor result in coupling it 
to an ergodic heat bath. In all six pairings of the equilibrium ergodic thermostats there is 
no evidence of dissipation or reduced phase-space dimensionality. Evidently all three (C;0 
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Figure 5; The eight equilibrium Lyapunov exponents describing pairs of ergodic thermostats 
linked by a Hooke’s-Law spring exhibit pairing when time averaged. The JB thermostats exhibit 
slightly larger Lyapunov exponents due to the quartic terms in their equations of motion. The 
summed spectra shown in the six upper lines give the growth rates of 1-, 2-, 3-, ... 8-dimensional 
phase volumes. The last sum is precisely zero as the equilibrium equations of motion conserve 
phase volume when time-averaged. 

thermostats are good heat baths from the standpoint of the Zeroth Law of Thermodynamics. 
Coupling pairs of thermostated oscillators with a harmonic coupling reveals that all three 
of the ergodic models { HH, JB, MKT } remain ergodic when coupled and also provide the 
complete canonical distribution. 

VI. NONEQUILIBRIUM SIMULATIONS 

The similar equilibrium behavior of the HH, JB, and MKT thermostats suggests a further 
comparison for simple heat-flow problems in which two coupled thermostats have different 
temperatures, leading to hot-to-cold heat flow and dissipation. Dissipation is reflected in 
phase space by the formation of strange attractors with a fractional dimensionality reduced 
below the equilibrium value. In the linear-response regime the six two-bath possibilities must 
necessarily agree with Green and Kubo’s theory as all the two-constant heat baths reproduce 
Gibbs’ canonical (g,p) distribution at equilibrium. Let us be adventurous and treat instead 
the relatively far-from-equilibrium coupling of two thermostats with thermostat 1 at half 
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Figure 6: The eight exponents describing the coupling of ergodic thermostats with markedly 
different temperatures, 1 and 2. Spectra for all nine combinations of thermostats are shown. The 
hot-to-cold heat current dissipates phase volume. On the average the phase volume approaches 
zero ( a strange attractor ) as exp[ X) despite the time-reversible nature of all the 

equations of motion. 

the temperature of thermostat 2 : 

pi = pi(T = 1) + ^2 - gi ; p2= p2{T = 2) + qi-q2 . 

Just as before all the masses, force constants, and relaxation times are chosen equal to unity. 

The far-from-equilibrium Lyapunov data corresponding to the nine combinations of ther¬ 
mostats, with temperatures of 1 and 2, are shown in Figure 6 , both individually and as 
sums. The broad maximum corresponding to sums of three or four exponents indicates 
that the maximum growth rate in phase space applies to four-dimensional volumes. The 
summed-up spectrum crosses zero near an abscissa value of 7 indicating that the dimen¬ 
sionality loss for this far-from-equilibrium problem is on the order of 1. Considerably larger 
losses, comparable to the number of particles, occur in “0^-model” chain simulations where 
each of a few dozen particles is tethered to its lattice site with a cubic restoring forced. 

Closer to equilibrium, with cold and hot temperature of 1.0 and 1.1, all nine thermostat 
combinations have a Lyapunov sum of order -0.01 while temperatures of 1.0 and 2.0 provide 
dissipative sums in the range from -0.20 to -0.30 . Compare Figures 6 and 7. Just as in 
heat conduction according to Fourier’s Law, we would expect a linear-response dissipation 
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Figure 7; The eight exponents describing the coupling of ergodic thermostats with similar tem¬ 
peratures, 1 and 1.1. All nine combinations are shown. The hot-to-cold heat current could be 
estimated with linear response theory and is essentially the same for all three thermostat types 
because the equilibrium phase-space distribution is the same Gaussian for all three thermostat 
types, f{q,p) oc . 

quadratic in the temperature difference driving the flow, varying as (Thot — Tcoiay ■ 

Loss of phase-space volume corresponds to thermodynamic dissipation ( through Gibbs’ 
relationship of phase volume to entropy, S = klnQ, where G is phase volume or number 
of states ). Coupling two baths, with the cooler bath first, gives the following “entropy 
production rates” : 

HH + HH ; 0.322 ; HH + JB : 0.254 ; HH + MKT : 0.213 ; 

JB + HH ; 0.299 ; JB + JB : 0.253 ; JB + MKT ; 0.207 ; 

MKT + HH ; 0.291 ; MKT + JB : 0.290 ; MKT + MKT ; 0.229 . 

We used quote marks above to remind the reader that simply summing the spectrum is not 
the same as averaging the dissipation for the cold and hot thermostats. 

Dissipation can also be reckoned in terms of the phase-space dimensionality of the dissipa¬ 
tive strange attractor, using the Kaplan-Yorke linear interpolation between the last positive 
dimension and the hrst negative one. For the nine heat-bath pairings the dimensionalities 
are : 

HH + HH ; 7.19 ; HH + JB : 7.32 ; HH + MKT : 7.38 ; 
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Figure 8: Coupling a Mexican Hat potential with twin well depths of 0.25 at g = ±1 to any of 
the three ergodic heat baths provides a mechanism for jumping over the barrier. The simulations 
shown here indicate that all three baths provide essentially the same temperature-dependent jump 
frequency v . The dependence is nearly linear when the coupling is small ( k = 1/10 ) and the 
abscissa is Gibbs’: . Eyring’s description ( the upper three curves ) is shown where the 

abscissa is now . Gibbs’ better approximates an Arrhenius straight line. Monte Carlo 

simulation of Mexican Hat dynamics, with a maximum jump length of 0.009\/T , results in the 
heavier line plotted in addition to the results for the three heat-bath models. 

JB + HH : 7.16 ; JB + JB : 7.25 ; JB + MKT : 7.33 ; 

MKT + HH : 7.02 ; MKT + JB : 7.01 ; MKT + MKT : 7.17 . 

These results should provide good benchmarks for nonlinear theories of transport. 


VII. JUMPS - EQUILIBRIUM MEXICAN HAT POTENTIAL SIMULATIONS 


The double-well “Mexican Hat” potential^ , 


0(g) = 


2 


+ 


(£_ 

4 


0(±1) 



0(0) = 0 , 


has a barrier between its two minima of height (1/4) . At a temperature of unity, where 
the barrier is negligible, simulations of the Hat potential coupled to the HH, MKT, and JB 
thermostats result in excellent agreement with Gibbs’ canonical distribution. At sufficiently 
low temperatures we expect that the frequency of successful jumps over the barrier will 
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decline, as would also be the case using the Metropolis, Rosenbluths, and Tellers’ Monte 
Carlo scheme with a reasonable jump length, perhaps of order 0.01 . 

A straightforward “weak coupling” of the two-minimum Hat potential to any one of the 
ergodic oscillator heat baths , 

‘h((?Bath, 9Hat) = —( Q'Bath “ Q'Hat Y i 

leads to “jump” frequencies that vary roughly as Gibbs’ canonical barrier weight, . 

Eyring’s model, as mentioned in the Wikipedia article on “Henry Eyring ( chemist )”, and 
with a jump frequency varying as , is a relatively poor description and corresponds 

to the upper curves in the Figure 8 . All three ergodic thermostat models behave similarly, 
as is shown in triplets of narrower lines in Figure 8 . 

Another possibility for modelling jumps is to use the Metropolis, Rosenbluths, and Tellers’ 
Monte Carlo method. By adusting the Monte Carlo trial steplength the two approaches, 
dynamical and Monte Carlo, can be made to correspond. We have compared the number of 
jumps over the barrier in billion-jump simulations with maximum jump length dq = 0.009\/T 
to dynamical simulations with timestep 0.005, where in all these latter dynamical cases the 
Hat potential is coupled to an HH, MKT, or JB thermostat with temperatures varying from 
0.010 to 1.00 , as shown in Figure 8 . 

VIII. SUMMARY 

The Lyapunov spectra for three varieties of time-reversible deterministic oscillator ther¬ 
mostats show that all three of them reproduce an extended ( four-variable ) form of Gibbs’ 
canonical ( two-variable ) distribution , 

The HH and MKT equations have Gaussian distributions for ( and ^ while the JB distri¬ 
bution is Gaussian in and All three approaches are sufficiently robust to thermostat a 
harmonic oscillator. The three are also consistent with the Zeroth Law of Thermodynamics 
in the sense that any pair of them at a common temperature T generates Gibbs’ canonical 
distribution for that temperature for both thermostats. 

This set of ergodic thermostats represents a good toolkit for simulations with a few 
degrees of freedom as it gives an idea of the dependence of dissipation on the chosen form of 
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the thermostat. With many degrees of freedom, where ergodicity is not an issue [ Poincare 
recurrence takes forever so that the size of fluctuations determines the usefulness of the 
results ] any of these thermostats, as well as the single thermostat Nose-Hoover model, are 
likely equally useful. 

We have considered two nonergodic time-reversible thermostats, the Patra-Bhattacharya 
and Nose-Hoover, which reproduce the specihed kinetic temperature but do not reproduce 
all of Gibbs’ canonical distribution. When the PB or NH thermostats are coupled to the HH, 
JB, or MKT thermostats two qualitatively different results occur: [ 1 ] a dissipative strange 
attractor in the case of the PB model and [ 2 ] a conservative integer-dimensional chaotic sea 
when coupled to the Hamiltonian-based NosGHoover oscillator. These results are similar 
for each of the three ergodic thermostats. In the PB case the dissipation is an indicator 
of the nonHamiltonian nature of the dynamics. Any of the ergodic thermostats, Hoover- 
Holian, Ju-Bulgac, Martyna-Klein-Tuckerman, can be used successfully in both equilibrium 
and nonequilibrium simulations. 

The hrst two of these can also be applied in nonequilibrium simulations where it is 
desired to specify the kinetic temperature of selected degrees of freedom. Though there is 
no problem in solving the MKT equations for nonequilibrium problems Brad Holian pointed 
out long agoi^ that this thermostat fails to return its target temperature due to the nonzero 
correlation linking the two MKT friction coefficients, ( C'C )mk^ 7^ 0 • 

This relatively thorough investigation of ergodicity, based on one million independent 
initial conditions, should settle ( at least from the numerical, as opposed to analytical, 
viewpoint ) the question of the ergodicity of the Martyna-Klein-Tuckerman oscillator. The 
apparent “holes” in that oscillator’s cross section are still a small “puzzle”. A second such 
“puzzle” is the nature of the local Lyapunov exponents. Their intricate multifractal Lya¬ 
punov structure is well concealed within an innocent Gaussian distribution. 

From the standpoint of dynamical-systems theory ( as opposed to thermodynamics ) 
we wish to emphasize the amazing nature of the chaos hidden within the simple ergodic 
Gaussian distributions. The smoothness of the distributions and the good convergence of 
their moments conceals the fractal nature of their chaos, as shown in Figure 1 through 3 . 
We recommend all three systems for further studies. The (C = ^ = 0) plane is only a single 
choice among the many possible. Despite their fractal nature the cross-sectional values are 
well-behaved, accessible, and reproducible from the numerical standpoint. It is encouraging 
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for the future that work originally chosen to settle some ergodicity questions turned out to 
generate a swarm of other problems still needing more detailed understanding. 
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